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We introduce a simple method for characterizing reactive pathways in quantum systems. Flux auto- 
correlation and cross-correlation functions are employed to develop a quantitative measure of dynamical 
coupling in quantum transition events, such as reactive tunneling and resonant energy transfer. We utilize 
the method to study condensed-phase proton-coupled electron transfer (PCET) reactions and to determine 
the relative importance of competing concerted and sequential reaction pathways. Results presented here 
include numerically exact quantum dynamics simulations for model condensed-phase PCET reactions. This 
work demonstrates the applicability of the new method for the analysis of both approximate and exact 
quantum dynamics simulations. 



I. INTRODUCTION 

Detailed mechanistic understanding of condensed- 
phase chemical dynamics is essential for the design of 
next-generation molecular catalysts and photosystems. 
The role of numerical simulations in this effort is partic- 
ularly important for systems that exhibit multiple, com- 
peting reactive pathways. The characterization of reac- 
tive pathways in classical systems has been greatly ad- 
vanced by the development of methods for rare event 
samplingl"2l. However, the corresponding tools for quan- 
tum systemjpHIIl are less developed, despite significant 
demand for understanding reactive charge and energy 
transfer pathways in complex quantum systems. 

In this paper, we introduce a method for characteriz- 
ing reactive pathways and quantifying dynamical corre- 
lations in general quantum systems using real-time flux 
auto-correlation (FAC) and flux cross-correlation (FCC) 
functions. We employ the new approach to investi- 
gate the reaction dynamics of condensed-phase proton- 
coupled electron transfer (PCET), an important proto- 
type for quantum systems that exhibit multiple reaction 
pathwaysSSHIIl Numerical results demonstrate the prac- 
ticality and utility of the new method for the mechanistic 
analysis of coupled quantum dynamical processes. 



II. THEORY 

We consider a set of states {\n)} that span the full 
Hilbert space of a closed quantum system, and we con- 
sider the partitioning of this set of states into N subsets 
{Qj}, such that each state belongs to one and only one 
subset. The criteria for this partitioning are flexible; for 
instance, cut-offs based on energy expectation values can 
be used to generate subsets of states with similar ener- 
gies, and conditions based on position expectation values 
can be used to generate subsets of states with similar con- 
figurations. The dividing surface corresponding to subset 
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The net flux associated with transitions from subset 
ilj to the remainder of the Hilbert space is given by 

F J = l l [H,V ] ] 1 (3) 

where H is the Hamiltonian for the full system. Dynami- 
cal correlations between these transitions can be obtai ned 
from thermal, real-time FAC and FCC functions^HH, 

C mj (t) = Tr [e-P H F m e iHt F je - im ] , (4) 

where (3 is the reciprocal temperature, and the FAC func- 
tion corresponds to the case m — j. Physically, Eq. Q 
describes the correlation between transitions into/out-of 
subset f2 m with transitions into/out-of subset Qj. 

Using that the net flux for a closed quantum system 
is zero, i.e. 2j=i Fj = 0' tne FAC function can be ex- 
pressed as a linear combination of FCC functions; specif- 
ically, 

N 

C J3 (t) = -Y^C mi (t), (5) 

m— 1 

where the summation excludes the m = j term. Intro- 
ducing the zeroth moments of the real part of these cor- 
relation functions, 

k mj = / dtRe[C mj (t)}, (6) 
Jo 
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Eq. ([5]) then yields 



JV 



(7) 



Within the assumption that Vj divides the system into 
two basins of stability, Eq. Q provides a decomposition 
of the corresponding reaction rate. We thus introduce 
a measure of the degree to which dynamical correlations 
between transitions associated with subsets £lj and Q m 
contribute to the overall reaction rate, 



h :jj 



(8) 



More generally, regardless of metastability of subsets, the 
dynamical correlation factor (DCF) defined in Eq. ^ 
provides a transferable measure of the relative contribu- 
tion from subset fl m to the transient dynamics associ- 
ated with entering/leaving subset Clj. Non-zero values 
for the DCF indicate correlated transitions between the 
subsets and provide a basis for identifying important dy- 
namical pathways and reaction mechanisms. We note 
that the DCF can be calculated using either exact or 
approximate quantum dynamical methods^, and flexi- 
bility in the definition of subset partitions enables the 
detailed characterization of complex systems. Further- 
more, to characterize processes that are far from thermal 
equilibrium, DCFs can be similarly constructed in terms 
of flux correlation functions with non-Boltzmann initial 
distributions. 



III. PCET MODEL SYSTEMS 

To demonstrate this flux-correlation approach for the 
analysis of systems with competing dynamical pathways, 
we consider a series of condensed-phase PCET reactions. 
The reactions are described using a system-bath Hamil- 
tonian, 
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+ + V v (x) + V ps (x, s) + V e {x, s) + 
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where the coordinates x, s and Qj correspond to the pro- 
ton, solvent polarization, and bath modes with masses of 
m x , m s , and M~, respectively. Potential energy surfaces 
associated with the donor and acceptor electronic states 
are described in the diabatic representation, 

V(r h)- f Vll ^ + Vc ^ x) Vl2{s) \ (W) 

Ve[X ' S) ~{ V 12 (s) V 22 (s) - V cp (x) ) ' ( W > 

where Vu(s) = ^m s uj 2 {s — si) 2 with 
i = 1,2, the electron-proton coupling is given by 
V op (x) = (1,2 tanh(0x), and V\2{s) is the diabatic cou- 
pling. 



The term V p (x) in Eq. (l9| describes a double- well po- 
tential for the proton coordinate, 
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where uj x is the proton vibrational frequency, and 
Vq is the intrinsic proton transfer barrier height. 
The proton solvent interaction is described using 
V ps (x, s) = — fj,is tanh(0x). 

The dissipative bath in Eq. Q exhibits an Ohmic 
spectral density, 



(12) 



where uj c is the cut-off frequency, and 77 is the dimension- 
less coupling between solvent and bath moded^. The 
spectral density is discretized using / oscillators with 
frequencies^ 



LJj = -uj c log 



3 - 0-5 
/ 



and coupling constants 



f 2riM i oj c 



(13) 



(14) 



where j = 1, . . . /. 

PCET reactions proceed via concerted or sequential 
mechanisms, depending on the chronology of the electron 
transfer (ET) and proton transfer (PT) events. The con- 
certed mechanism involves dynamically correlated trans- 
fer of both the electron and proton, whereas the sequen- 
tial mechanism involves dynamically uncorrelated ET 
and PT steps. We consider three models for PCET in 
the current study, each of which corresponds to a dif- 
ferent set of parameters for the Hamiltonian in Eq. J9J. 
Models I and II are adapted from an earlier study of uni- 
directional PCET in non-dissipative systems and assume 
that the solvent responds to the transferring electron and 
proton as if they are like-charged particles^. This un- 
physical picture of solvent polarization is corrected in 
model III, in which the solvent response to the trans- 
ferring electron is counteracted by the solvent response 
to the transferring proton. The parameters for all three 
models are provided in Table [i] 

IV. CALCULATION DETAILS 

For each PCET model system, FAC and FCC functions 
are calculated usingthe quasi-adiabatic path-integral 
(QUAPI) methocPlEl]. QUAPI is a numerically exact 
quantum dynamics method that employs a real-time path 
integral formulation. It has previously been used to study 
single-particle transfer reactions, such as ET or 
here, we extend the approach to describe the coupled 
transfer of both an electron and a proton. Details of the 
QUAPI implementation are provided in Appendix I. 
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TABLE I. Parameters for the model PCET systems. 



Parameter a Model I Model II Model III 





1836.1 


1836.1 


1836.1 




0.0104 


0.0104 


0.0104 


V 


0.012 


0.014 


0.012 


m 3 


22000 


22000 


22000 


LJ S x 10 4 


3.72 


4.00 


3.72 


Si 


-2.13 


-2.16 


-2.13 


S2 


2.13 


2.16 


2.13 


Vl2 


0.0245 


0.0124 


0.0245 


A*i 


0.0011 


0.017 


-0.0011 


/i 2 x 10 3 


5.84 


0.71 


5.84 


<t> 


8.0 


8.0 


8.0 


A 


0.0 


0.012 


0.0 


/ 


12 


12 


12 


M :l 


m a 


rn s 


m s 


V 


m s cu s 


m s uj s 


m s u) s 


T/K 


300 


750 


300 



a All values in atomic units, unless otherwise specified. 

We partition the full set of quantum states for the 
PCET Hamiltonian (Eq. ^) into N = 3 subsets that are 
defined in terms of the proton position and the electronic 
diabatic state. The subset associated with the PCET re- 
actant states, Sl r , is defined in terms of the projection 
operator 

V t = \l){l\h(-x). (15) 

Similarly, f2 p includes the PCET product states and is 
defined using 

V p = \2)(2\h(x), (16) 

and f2; includes intermediate states associated with 
single-particle transfer and is defined using 

A = |l><l|/i(ar) + |2><2|M-a:). (17) 

In these equations, |1) and |2) indicate the donor and 
acceptor electronic diabatic states, and h(x) is the heav- 
iside function 

Using these subset definitions, the corresponding FAC 
and FCC functions are calculated using Eq. Q. Finally, 
Eq. |8| is used to calculate k\ p and K rp , which quantify 
dynamical correlations of the product subset with the 
intermediate and reactant subsets, respectively. For the 
present case in which the states are simply partitioned 
into reactant, product, and intermediate subsets, we can 
relate the calculated DCF to competing reaction mech- 
anisms. Specifically, by reporting on whether the reac- 
tion dynamics proceeds via the intermediate subset, or 
whether it proceeds via direct transfer from the reactant 
to product subsets, K- lp and K rp respectively indicate the 
degree to which the sequential or concerted PCET mech- 
anism is dominant. 



V. RESULTS 

Fig. [TJ presents the calculated FAC and FCC functions 
for model I, with C pp (t), C rp (i), and C\_ p (t) plotted in 
blue, red, and green, respectively. The FAC function ex- 
hibits initial decay on the 500 a.u. timescale, followed by 
modest recrossing. Cross-correlations in the subset dy- 
namics are most pronounced for the reactant and product 
subsets, with C rp (t) closely mirroring the features of the 
product FAC function. In contrast, only a small degree 
of cross-correlation is found for the dynamics associated 
with the intermediate and product subsets. All dynam- 
ical correlations among the subsets are found to vanish 
by approximately 4000 a.u. in time. 




time (a.u.) 



FIG. 1. Numerically exact, symmetrized flux-correlation 
functions for model I, with C pp (i) in blue, C rp (t) in red, and 
Cip(t) in green. 

Integration of these correlation functions yields the 
DCF for model I, which are reported in Table [Tl] The 
larger magnitude of K rp indicates that the PCET reaction 
in model I primarily proceeds via the concerted mecha- 
nism, although a substantial contribution from the se- 
quential pathway is also found. We note that earlier 
simulations for a non-dissipative version of this PCET 
model als o con cluded the importance of the concerted 
mechanis m 1 26 1 32 1 . 

TABLE II. Reaction mechanisms for PCET in models I-III. 
^Pathway DCF Model I Model II Model TIT 
Concerted K Tp 0.62(6) 0.010(7) 0.93(7) 
Sequential m p 0.38(5) 0.99(5) 0.07(3) 



Fig. [2] presents the C pp (t), C rp (t), and Ci p (t) correla- 
tion functions for model II at short times, which reveal 
significant differences in comparison to Fig. [T] The FAC 
function for model II exhibits pronounced oscillations on 
the timescale of proton vibrations, but more striking dif- 
ferences are seen in the cross-correlation functions. In 
Fig. [21 far greater contributions are seen from the C lp 
than o rp , indicating that flux into the product state is 
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dynamically coupled with the single-particle transfer in- 
termediates, rather than with the reactant subset. This 
point is further illustrated in Table |TlJ which reveals «j p 
to approach unity while K rp nearly vanishes. These re- 
sults clearly indicate that model II exhibits a sequential 
PCET reaction mechanism, in a greem ent with studies of 
a related, non-dissipative systenpfil22l. 




time (a.u.) 



FIG. 2. Numerically exact, symmetrized flux-correlation 
functions for model II, with C pp (t) in blue, C Tp (t) in red, and 
Cip(i) in green. 

Finally, Fig. [3] presents the flux-correlation functions 
for model III, which differs from model I only in terms 
of the solvent polarization response to the proton coordi- 
nate (Table [TJ) . The FAC function for model III exhibits 
a slower timescale for initial decay than for model I, as 
well as a more pronounced degree of dynamical recross- 
ing. Both FCC functions reflect the time dependence 
of the FAC function, although the magnitude of C rp (i) 
is greater than that of C- lv at all times. Interestingly, 
Table [Tl] shows that K lp approaches unity, whereas k; p 
nearly vanishes, an indication that the PCET reaction 
in model III is more dominated by the concerted mecha- 
nism than the corresponding concerted reaction in model 
I. In model III, the solvent response to the net-neutral 
PCET charge-transfer reaction yields reactant and prod- 
uct states with less solvent polarization than in model 
I; the energetic favorability of minimizing solvent reor- 
ganization throughout the reaction thus creates a driv- 
ing force for co-localizing the charge distributions of the 
transferring electron and proton, which leads to the more 
strongly concerted PCET mechanism observed for model 
III. 




i i i i i i i i i 

1000 2000 3000 4000 

time (a.u.) 



FIG. 3. Numerically exact, symmetrized flux-correlation 
functions for model III, with C pp (i) in blue, C rp (i) in red, 
and Cip(t) in green. 



phase PCET. In particular, we quantify the relative im- 
portance of the sequential and concerted PCET reaction 
mechanisms for two previously studied model systems, 
and in a third system, we demonstrate that introduction 
of a refined description for the solvent-proton coupling 
leads to greater dominance of the concerted PCET mech- 
anism. It is expected that this approach will prove use- 
ful in the future analysis of reaction and energy transfer 
pathways in complex systems. 
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APPENDIX I: QUAPI IMPLEMENTATION FOR PCET 



VI. CONCLUSIONS 

We introduce a method in which flux-correlation func- 
tions are used to characterize dynamical correlations and 
reaction pathways in quantum systems. Numerical re- 
sults demonstrate the utility of the method for the anal- 
ysis of exact quantum dynamic simulations of condensed- 



Here, we describe the implementation of the QUAPI 
method used in the current study. The total PCET 
Hamiltonian in Eq. ^ is written as a sum of system 
and bath contributions such that, H = H s + Hg, where 

2 2 

H s = ^ + ^+ V e {x, s) + V p (x) + V ps (x, s) (19) 
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and the Hamiltonian describing the bath modes and their 
coupling to the solvent coordinate is 



/ p2 f ^ 



2Mi 
=1 J 



t 2 M ri\?*-Mj3j) ' (20) 



In practice, we report the symmetrized form of the flux- 
flux correlation functions and note that the zeroth mo- 
ment of these functions is identical to that of the stan- 
dard correlation function defined in Eq. Q 21 . The sym- 
metrized FCC function is given by 



C im (t) = Tr^ 



JHfJh 



-iHt 



(21) 



where the subset indices j, to € {i,r,p}, and 
t c = t — i/3H/2. The complex-time propagators in Eq. 



(21 ) are discretized into Af time slices of length Ai c , and 
the trace over bath modes is evaluated analytically to 
yield 



2 

ds dx E] ^( s ) K(s, x, <x; t c ) 

"* J tr=l 



X (sx,Xl,a 1 \Fj\s2M+2,XZbr+2,<?ZSr+2) (22) 
X (s/V+2, ZjV+2, VAf+2\F m \sM +1 , Xtf+l, o-jv+i), 

where s = {si, . . . , S2M+2}, x = {xi, . . . ,x 2 Af+2}, and 
cr = {ax, . . . ,0-2^+2}- Numerical evaluation of the in- 



tegrals in Eq. ([22j is performed using two independent 
path-integral Monte Carlo simulations (MCI and MC2), 
as we have previously described in detaiP^. 

The term K (s, x, <r; i c ) in Eq. [22] is the path-integral 
representation for the complex-time propagators of the 
system Hamiltonian, 

M+l 

K(s,x,a-;t c ) = 

k=2 

2A/+2 

x Yl (sk,x k ,o-k\e iHBM *\s k _ 1 ,x k _ 1 ,a k _ 1 },(23) 

k=Af+3 

where 



Y[ (s k ,x k ,<r k \e lH3 ^ tc \s k -i,x k -i,cj k -i) 



(s k ,x k ,a k \e iHsAi< =/ ft |s fc _i,a; fc _i,o- fe _i) 



(24) 



M 

E 

m— 1 



0m(sfc,^fc,c r fe)0^(s fc _i,a: fc _i,cr fe _i)e jE ™ A W^ 



4> m and E m are the eigenstates and eigenenergies of 
H s , respectively, and Mq is the number of eigen- 
states included in the expansion. The eigenstates 
and eigenvalues are obtained from a three-dimensional 
DVR grid calculation in terms of the solvent coordinate 
{— s* , — s* + As, . . . , s* — As, s*}, the proton coordinate 
{— x* , —x* + Ax, . . . , x* — Ax, x*}, and the two electronic 
states. 

The discretized form of the non-local influence func- 
tional in Eq. [22] is 

(2M+2 fe \ 
-EE B kk'S k Sfc' J , (25) 
fc=l fe'=l / 



where Iq is the p artition function of the uncoupled bath 
osciUator^ZlMiasi. For the case of linear system-bath cou- 
pling, the diagonal matrix elements are given by^ 



B kk = e 



/ c 2 



uj{t k+ i -t k + i(3) 



sin 



Uj(tk+1 — t k ) 



x sin 



(26) 



B kk '= E 



and the off-diagonal matrix elements are given by 

' UJj{t k +l — t k ) 

x cos 



- tk'+i + t k - t k > + iP) 



Uj(t k '+i — t k i) 



(27) 



The complex times t k in Eqs. (26) and (27) are provided 
in Table Hm 

TABLE III. Complex times t k used to calculate the {B kk /}. 



k 




t k 


1 







2,...,JV+1 


(*- 


l/2)At 


M + 2 


t - 


iph/2 


N +3,...,2N + 2 (2N - 


f 3/2 


- k)At* - iph 


27V + 3 




-iph 



The PCET flux operators, which correspond to the 
subset projection operators defined in Eqs. (15)-(17), are 



F p = F x \2)(2\ + h(x)F c , (28) 
F x = -Fj(l){\\-h{-x)F Bt (29) 
Fi = F x (|1)(1| - |2)(2|) + (h(-x) - h(x)) F c , (30) 

where F^ = i [H, h(x)] is the flux operator for the proton 
coordinate and F D = ^[H,\2)(2\] is the flux operator 
for the electronic diabatic states. The matrix elements 



appearing in Eq. ([28 ) are evaluated using 

(2\s k >,x k ',a k >) = 

Sa k ,2 $<tu,2 S(s k - S k >) 



(sfc , x k , at I F x 1 2) (2 1 s k > , x k > , crfc' ) = 
ill 



2m x xfd 

x [5(x k )S(x k f - x FB ) - 5(x k - x F D)S{x k >)} , 
where xfd = Ax, and 

(s fe , x k , a k \h(x)F e \sk' , x k > , <j k <) = 



(31) 



-h{x k )8{x k - x k ,)V 12 {s k )6(s k - s fe ') 



X <W Sa„ 



(32) 



respectively. Similar expressions are used for the terms 
in Eqs. H and g. 



All convergence parameters for the QUAPI calcula- 



tions are provided in Table IV Convergence with respect 
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to the number of path-integral beads is determined by 
comparing FAC functions with J\f — 3 — 5; the employed 
value of Af = 4 is consistent with the number for beads 
required in previous QU API simu lations for proton and 
electron transfer reaction d 29 * 31 * 35 l Convergence with re- 
spect to the eigenfunction expansion is determined by 
comparing the trace over the complex-time propagator 
with Mq = 100 — 2000; no significant numerical changes 
are found for cutoffs larger than those reported in Table 
IVl 

TABLE IV. Convergence parameters for QUAPI. 



Model 


Mo 


MCl a 


MC2 a 


s* 


As 


.X'* 


Aa; 


I 


500 


1 x 10 7 


1 x 10 s 


5 


0.1 


4 


0.2 


II 


750 


1 x 10 8 


5 x 10 s 


11 


0.09 


1 


0.14 


III 


500 


1 x 10 7 


1 x 10 s 


5 


0.1 


4 


0.2 



a Indicates the number of Metropolis acceptance/rejection steps. 

For numerical comparison, the rate of the concerted 
PCET reaction in model I is calculated using both the 
QUAPI method and a golden-rule rate expression. The 
QUAPI rate is obtained using 



k 



QUAPI 



"TP 

Qt' 



(33) 



where fc rp is defined in Eq. ^ and Q r is the parti- 
tion function for reactant subset f2 r ; this yields a rate 
feq = 1.7(1) x 10~ 6 a.u. The golden-rule approxima- 
tion for vibronically nonadiabatic PCET is given by^^ 



«GR 



V 2 



xexp 



7r\ 

■P(AG° + A + e„ 



4A 



. (34) 



where e M and e„ are the energies of the proton vibra- 
tional states \i and u, respectively, and P M is the ther- 
mal probability corresponding vibrational state /i. The 
PCET reaction in model I is ground-state dominated, 
electronically adiabatic and vibrationally nonadiabatic; 
construction of the diabats for model I thus yields a vi- 
bronic coupling of V 00 = 3.75 x 10~ 4 , solvent reorgani- 
zation energy A = 1.34 x 10 -2 a.u., and driving force 
AG = 0. The golden-rule rate obtained for the con- 
certed reaction in model I is /cqr = 2.06 x 10~ 6 a.u., 
which is in good agreement with the QUAPI rate stated 
above. 
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